// #include "opencv2/core.hpp"
#include <opencv2/opencv.hpp>
#include <opencv2/ximgproc.hpp>
#include <iostream>
#include <time.h>

using namespace cv;
using namespace std;

void gabor_filter(bool if_opencv_kernal,Mat gray_img, Mat& gabor_img, Mat gabor_tmp,int k=9, float sigma=1.0,float gamma=0.5,float lambda=5.0,float psi=-CV_PI / 2);
Mat gabor_kernal_wiki(Size ksize,double theta,double sigma=1.0,
                            double lambd=5.0, double gamma=0.5, double psi=-CV_PI / 2, int ktype=CV_32F);

int main()
{   
   
    clock_t total_start=clock();
    Mat src=imread("1.jpg");
    Mat src_gray,src_gray_gabor,gabor_tmp;
    cvtColor(src,src_gray,COLOR_BGR2GRAY);

    bool if_opencv_kernal=true;//false;
    gabor_filter(if_opencv_kernal,src_gray, src_gray_gabor, gabor_tmp);
    
    char save_name[100];
    sprintf(save_name, "gabor_filter_img.jpg"); 
    imwrite(save_name,src_gray_gabor);
    
    clock_t total_end=clock();
    cout<<"Total Time:"<<(double)(total_end - total_start)/ CLOCKS_PER_SEC*1000 <<endl;
    return 0;

}


/*Gabor 滤波处理相关
Sigma : σ 表示高斯函数的标准差
Gamma： γ 表示长宽比，决定这Gabor核函数图像的椭圆率
Lambda：λ 表示滤波的波长
Psi  ： ψ 相位偏移量，取值范围是-180~180
Theta： θ  表示Gabor核函数图像的倾斜角度 
*/
void  gabor_filter(bool if_opencv_kernal,Mat gray_img, Mat& gabor_img, Mat gabor_tmp,int k, float sigma,float gamma,float lambda, float psi) 
{   
    int ddepth=CV_8U;  
    //表示Gabor核函数图像的倾斜角度 // 0  45 90 135度 的核
    double theta[4]={
        0.0,
        CV_PI/4,
        CV_PI/2,
        CV_PI/4*3,
    }; 

    Size ksize=Size(k,k); //卷积核的大小
    Mat gabor_kernel;
    char save_name[100];
    gabor_img=Mat::zeros(gray_img.size(),CV_8UC1);
    // 分别取 0  45 90 135 四个方向gabor核进行计算
    for  (int i=0;i<4;i++)
    {       
        cout<<"theta θ= "<<i*45<<"="<<theta[i]<<endl;
        if(if_opencv_kernal) 
            gabor_kernel= getGaborKernel(ksize, sigma, theta[i], lambda, gamma, psi, CV_32F);
        else 
            gabor_kernel=gabor_kernal_wiki(ksize,theta[i]);
        // cout<<gabor_kernel<<endl;
        filter2D(gray_img,gabor_tmp,ddepth,gabor_kernel);

        /*合成各个角度的图*/
        //cv::max(src1,src2,dst);
        cv::max(gabor_tmp,gabor_img,gabor_img);
        sprintf(save_name, "im_resize_gray_gabor_%d.jpg",i*45); 
        imwrite(save_name,gabor_img);
    }



}


Mat gabor_kernal_wiki(Size ksize, double theta,double sigma, 
                            double lambd, double gamma, double psi, int ktype)
{
    double sigma_x = sigma;   //  σ 表示高斯函数的标准差
    double sigma_y = sigma/gamma; // γ 表示长宽比，决定这Gabor核函数图像的椭圆率

    int xmin, xmax, ymin, ymax;
    // double c = cos(theta), s = sin(theta);

    xmax = ksize.width/2;
    ymax = ksize.height/2;
    xmin = -xmax;
    ymin = -ymax;

      //------这部分内容是为了确定卷积核的大小------

    CV_Assert( ktype == CV_32F || ktype == CV_64F );

    Mat kernel(ymax - ymin + 1, xmax - xmin + 1, ktype);
    // double scale = 1;   
    // double ex = -0.5/(sigma_x*sigma_x);
    // double ey = -0.5/(sigma_y*sigma_y);
    /*
    y
    xmin, ymin,xmax, ymax -4 -4 4 4
    [[-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]
    [-4 -3 -2 -1  0  1  2  3  4]] 
    x:
    [[-4 -4 -4 -4 -4 -4 -4 -4 -4]
    [-3 -3 -3 -3 -3 -3 -3 -3 -3]
    [-2 -2 -2 -2 -2 -2 -2 -2 -2]
    [-1 -1 -1 -1 -1 -1 -1 -1 -1]
    [ 0  0  0  0  0  0  0  0  0]
    [ 1  1  1  1  1  1  1  1  1]
    [ 2  2  2  2  2  2  2  2  2]
    [ 3  3  3  3  3  3  3  3  3]
    [ 4  4  4  4  4  4  4  4  4]]

    */
   for( int x = xmin; x <= xmax; x++ )
   {
    // cout<<x<<endl;
    for( int y = ymin; y <= ymax; y++ )
        {   
            // cout<<y<<" ";
            double x_alpha=x*cos(theta)+y*sin(theta);
            double y_alpha=-x*sin(theta)+y*cos(theta);
            // double xr = x*c + y*s;
            // double yr = -x*s + y*c;
            // double exponent=std::exp(-0.5*(x_alpha*x_alpha/sigma_x*sigma_x+
            double exponent=exp(-0.5*(x_alpha*x_alpha/pow(sigma_x,2)+
                                           y_alpha*y_alpha/pow(sigma_y,2)
                                     )
                                );
            double v = exponent*cos(2*CV_PI/lambd*x_alpha + psi);
            if( ktype == CV_32F )
                {
                // printf("%d %d %d %d %d %d ",ymax,y, x,xmax, x+xmax,y+ymax);
                // printf("%d %d %d %d ", ymax - y, xmax - x,x+xmax,y+ymax);   
                kernel.at<float>(y+ymax,x+xmax) = (float)v;  // 从
                // kernel.at<float>(ymax - y, xmax - x) = (float)v;  // opencv的赋值
                }
            else
                kernel.at<double>(y+ymax,x+xmax) = v;
        }
    // cout<<endl;
   }

    cout<<kernel<<endl;
    return kernel;
}
